function out = yieldCurveDecom(setup,output,theta2)
% Information: we use the value of sigma for which theta11 is estimated to compute the bond yields.
% To simulate the factors, we use the value of sigma from the VAR model in step 3
% compute TP for all maturities
tmp                = output.model;
nm                 = setup.nm;
setupAll           = setup;
setupAll.matSelect = 1:max(setup.matSelect);
if output.model.modelNum == 1
    [modelAll,errorMes]= solveQTSM(tmp.alpha,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setupAll);
elseif output.model.modelNum == 2
    [modelAll,errorMes] = solveQTSM_constRsmooth(tmp.alpha,tmp.rhor,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setup);
end
modelAll.modelNum = output.model.modelNum;
modelAll.factorSignResOn = 0;
modelAll.macros   = output.macros;
modelAll.rLag     = output.model.rLag;
modelAll.rhor     = output.model.rhor;

% Getting the P dynamics. Here we use the estimates from the VAR model in step 3
factors          = [output.macros'; output.factors];
[nx,T]           = size(factors);
[h0,~,hx,~,sigmaVar] = unfoldTheta2_threshold_macro(theta2,nx);
sigmaVarSq       = sigmaVar*sigmaVar';
maxMat           = max(modelAll.matSelect);
%% Expected factor moments
% Expected factors under P
xExp = nan(nx,maxMat,T);
for t = 1:T
    for i = 1:maxMat
        if i == 1
            xExp(:,1,t) = factors(:,t);
        else
            xExp(:,i,t) = h0 + hx*xExp(:,i-1,t);
        end
    end
end

% Factor covariances under P
xCov = zeros(nx,nx,maxMat);
for i=1:maxMat
    if i == 1
        % xCov = 0 here
    elseif i == 2
        xCov(:,:,i) = sigmaVarSq;
    elseif i > 2
        xCov(:,:,i) = xCov(:,:,i-1) + hx^(i-1)*sigmaVarSq*(hx^(i-1))';
    end
end

%% Expected short rates
if modelAll.modelNum == 1
    % Expected short under P
    rExp = nan(maxMat,T);
    for t = 1:T
        for i = 1:maxMat
            rExp(i,t) = modelAll.r0 + modelAll.rx*xExp(:,i,t) ...
                + modelAll.rxx*kron3(xExp(:,i,t),xExp(:,i,t)) ...
                + trace(reshape(modelAll.rxx,nx,nx)*xCov(:,:,i));
            % alternatively
            %rExp(i,t) = modelAll.r0 + modelAll.rx*xExp(:,i,t) ...
            %    + modelAll.rxx(:)'*reshape((xExp(:,i,t)*xExp(:,i,t)'+xCov(:,:,i)),nx^2,1);
        end
    end
elseif modelAll.modelNum == 2
    rExp = nan(maxMat,T);
    for t = 1:T
        for i = 1:maxMat
            if i == 1
                rExp(i,t) = modelAll.r0 + modelAll.rhor*modelAll.rLag(t,1) ...
                    + modelAll.rx*xExp(:,i,t) + modelAll.rxx*kron3(xExp(:,i,t),xExp(:,i,t)) ...
                    + trace(reshape(modelAll.rxx,nx,nx)*xCov(:,:,i));
            else
                rExp(i,t) = modelAll.r0 + modelAll.rhor*rExp(i-1,t) ...
                    + modelAll.rx*xExp(:,i,t) + modelAll.rxx*kron3(xExp(:,i,t),xExp(:,i,t)) ...
                    + trace(reshape(modelAll.rxx,nx,nx)*xCov(:,:,i));
            end
        end
    end    
end

%% Fitted yields and term premia
% tp(n)_t = y(n)_t - 1/n*sum(E_t[r_t+i],i=0,1,...n-1)
ny         = size(modelAll.matSelect,2);
yHat       = nan(ny,T);
termPrem   = nan(ny,T);
rExpAvg    = nan(ny,T);
for t=1:T
    yHat(:,t) = modelFit(modelAll,factors(nm+1:end,t),t);
    for i=1:ny
        rExpAvg(i,t) = mean(rExp(1:modelAll.matSelect(1,i),t));
        termPrem(i,t) = yHat(i,t) - rExpAvg(i,t);
    end
end

%% Convert into forwards
fHat    = nan(ny,T);
fwdPrem = nan(ny,T);
for i = 1:ny-1
    fHat(i,:)    = (i+1)*yHat(i+1,:) - i*yHat(i,:);
    fwdPrem(i,:) = (i+1)*termPrem(i+1,:) - i*termPrem(i,:);
end

%% Output
out.yieldsHat     = yHat;
out.termPremia    = termPrem*10000;
out.rExp          = rExp;
out.forwardsHat   = fHat;
out.forwardPermia = fwdPrem*10000;
out.matSelect     = modelAll.matSelect;


end